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Abstract 

Global sensitivity analysis is often impracticable for complex and time demanding numerical mod- 
els, as it requires a large number of runs. The reduced-basis approach provides a way to replace 
the original model by a much faster to run code. In this paper, we are interested in the information 
loss induced by the approximation on the estimation of sensitivity indices. We present a method 
to provide a robust error assessment, hence enabling significant time savings without sacrifice on 
precision and rigourousness. We illustrate our method with an experiment where computation 
time is divided by a factor of nearly 6. We also give directions on tuning some of the parameters 
used in our estimation algorithms. 

Keywords: sensitivity analysis, reduced basis method, Sobol indices, bootstrap method, Monte 
Carlo method. 



Introduction 

Many mathematical models use a large number of poorly-known parameters as inputs. When such 
models are encountered, it is important for the practitioner to quantify whether this uncertainty on 
the inputs has a large repercussion on the model output. This problem can be tackled by turning 
the uncertain input parameters into random variables, whose probability distribution reflects the 
practitioner's belief about the oddness of the fact that an input parameter takes some value. 
In turn, model output, as function of the model inputs, is a random variable; its probability 
distribution, uniquely determined by the inputs' distribution and the model itself, can give detailed 
and valuable information about the behavior of the output when input parameters vary: range of 
attained values, mean value and dispersion about the mean (throughout expectation and standard 
deviation), most probable values (modes), etc. 

Sensitivity analysis aims to identify the sensitive parameters, that is the parameters for which a 
small variation implies a large variation of the model output. In stochastic sensitivity analysis, one 
makes use of the output's probability distribution to define (amongst other measures of sensitivity) 
sensitivity indices (also known as Sobol indices). Sensitivity index of an output with respect to 
an input variable is the fraction of the variance of the output which can be "explained" by the 
variation of the input variable, either alone (one then speaks about main effect), or in conjunction 
with other input variables (total effect). This way, input variables can be sorted by the order 
of importance they have on the output. One can also consider the part of variance caused by 
the variation of groups of two or more inputs, although main effects and total effects are generally 
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sufficient to produce a satisfying sensitivity analysis. The reader is referred to Helton et al. (2006); 
Saltelli et al. for more information about uncertainty and sensitivity analysis. 
Once these indices have been defined, the question of their effective calculation remains open. For 
most models, an exact, analytic computation is not attainable (even expressing an output as an 
analytic function of the inputs is infeasible) and one has to use numerical approximations. 
A robust, popular way to obtain such approximations is Monte Carlo estimation. This method 
simulates randomness in inputs by sampling a large number of parameters' values (from the selected 
inputs' distribution). The model output is then computed for each sampled value of the parameters. 
This way, one obtains a sample of outputs, under the conjugate action of the model and the input 
distribution. This sample of outputs is fed into a suitable statistical estimator of the desired 
sensitivity index to produce a numerical estimate. The Monte Carlo approach to computation 
of Sobol indices is described in Sobol (2001), together with improvements in Homma and Saltelli 
(1996); Saltelli (2002). 

A major drawback of the Monte Carlo estimation is that a large number of outputs of the model 
have to be evaluated for the resulting approximation of the sensitivity index to be accurate enough. 
In complex models, where a simulation for one single value of the parameters may take several 
minutes, the use of these methods "as-is" is impracticable. In those cases, one generally makes 
use of a surrogate model (also known as reduced model, metamodel or response surface). The 
surrogate model has to approximate well the original model (called the full model), while being 
many times faster to evaluate. The sensitivity index is then calculated via a sample of out- 
puts generated by a call to the surrogate model, thus accelerating the overall computation time. 
The reduced-basis (RB) method Nguyen et al. (2005); Grepl and Patera (2005); Veroy and Patera 
(2005); Grepl et al. (2007) is a way of defining surrogate models when the original model is a dis- 
cretization of a partial differential equation (PDE) depending on the input parameters. It comes 
with an error bound, that is, an upper bound on the error between the original output and the 
surrogate output. 

The sensitivity index produced by Monte Carlo estimation on a surrogate model is tainted with 
a twofold error. Firstly, our Monte-Carlo sampling procedure assimilates the whole (generally 
infinite) population of possible inputs with the finite, randomly chosen, sample; this produces 
sampling, or Monte-Carlo error. Secondly, using a surrogate model biases the estimation of the 
Sobol index, as what is actually estimated is sensitivity of surrogate output, and not the full one; 
we call this bias the metamodel error. 

In order to make a rigorous sensitivity analysis, it is important to assess the magnitude of these 
two combined errors on the estimated sensitivity indices. Such assessment can also be used to 
help in the choice of correct approximation parameters (Monte-Carlo sample size and metamodel 
fidelity) to achieve a desired precision in estimated indices. 

Sampling error can be classically estimated for a moderate cost by using bootstrap resampling 
Efron et al. (1993); Archer et al. (1997). Based on statistical estimation theory, the bootstrap 
technique involves the generation of a sample of replications of sensitivity index estimator, whose 
empirical distribution serves as approximation of the true (unknown) estimator distribution, in 
order to produce asymptotic confidence intervals which give good results in many practical cases. 
A variation on the bootstrap, which addresses sampling error as well as metamodel error, has 
been proposed in Storlie et al. (2009); also Marrel et al. (2009) develops a methodology in Kriging 
metamodels. In this paper, we present another confidence interval-based approach for assessing 
sampling errors, together with errors caused by reduced-basis metamodels, which makes use of the 
certified, a posteriori error bound that comes with the reduced-basis method. Boyaval et al. (2009) 
also makes use of the reduced-basis output error bound to certify computation of the expectation 
and the variance of a model output with neglected sampling error. 

The advantages of our approach are: its rigorousness (the impact of the use of a surrogate model 
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is provably bounded), its efficiency (our bounds are rather sharp, and go to zero when metamodel 
errors decrease), its clear separation between estimation (sampling) and metamodel error, and 
moderate computational requirements (time should rather be spent at making a precise compu- 
tation than at measuring precision). In other words, our method allows to estimate sensitivity 
indices by using a reduced basis metamodel which largely speeds up computation times, while 
rigorously keeping track of the precision of the estimation. 

This paper is organized as follows: in the first part, we go through the prerequisites for our 
approach: we give the definition and standard Monte Carlo estimator of the sensitivity indices we 
are interested in, and give an overview of the reduced basis method; in the second and third parts, 
we present our confidence interval estimation technique for the sensitivity index, which accounts 
for the two sources of error described earlier. In the fourth part, we present the numerical results 
we obtain on an example of a reduced-basis metamodel. 



1. Model output and sensitivity analysis methodology 

1.1. Sensitivity indices 

1.1.1. General setting 

In order to define sensitivity indices, we choose a probability distribution for the input variables, 
turning each input variable Xi {i = l,...,p) into a random variable with known distribution; 
the model output Y = f(Xi, . . . ,X p ) (assumed to be a scalar; multiple outputs can be treated 
separately) is thus for X = (X 1; . . . ,X p ) a cr(X)-measurable random variable. We further assume 
that the X^s are mutually independent. We also fix throughout all this paper an input variable 
of interest 1 < % < p. We define the first-order main effect of Xi on Y by: 

_ VarE(F|X t ) 
bi ~ VarF (lj 

Si is the sensitivity index in which we are interested in this paper but other indices (total effect, 
high-order effects) exist and our methodology can readily be extended to these indices. 

1.1.2. Monte-Carlo estimator 

We are interested in the following Monte-Carlo estimator for Si Homma and Saltelli (1996); Saltelli 
(2002): a sample size JVgN being given, let and |x' fc |^ ^ be two random i.i.d. 

samples of size N each, drawn from the distribution of the input vector X. 
For k = 1 N, we note: 



y k = /(x fc ) 

and: 

Vk — J l A l 5 • • • 5 A i-1> A j > A i + 1> • • • > A p ) 

The Monte-Carlo estimator of Si is then given by: 

-q~ _ ~k Sfc=i VkUk ~ {jj SfcLi Ukj yjj J2k=i y'k) 

N 2^k=l ilk \n ^k=l ilk 

It can be shown that Si is a strongly consistent estimator of Si. 
1.2. Metamodel construction: overview of the reduced basis method 

We briefly present here the reduced basis method for affinely parametrized elliptic partial differ- 
ential equations. For more details, see Nguyen et al. (2005). 
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1.2.1. Offline-online decomposition 

The reduced basis method deals with variational problems of the form: given an input parameter 
vector X belonging to a parameter set V C R p , 

find u(X) G 7 so that a(u, v; X) = ip(v) \/v G 7 (3) 

where T is an appropriate finite dimensional function space (usually a discretization of a continuous 
function space such as H 1 or Hq), i/j is a linear form on J 7 , and a(-, •; X) is a parameter-dependent 
bilinear form on J 7 that can be written under affine form: 

Q 

a(v, w; X) = 53 Q q (X.)a q (v, w) Vi>, w e J 7 

9=1 

where g are arbitrary real functions defined on X>, and a g are parameter-independent bilinear 
functions on J 7 . 

Traditional computation of u(X) for a prescribed X involves looking for u(X) as a linear combi- 
nation: 

dim T 

u(X) = 53 Mi (X)0, 
i=i 

where (0j) i=1 dimj 7 * s a suitable basis of J 7 , and the unknowns are (iti(X))j = i ) ... j dii n .F- This way 
one obtains the following linear system of (dim J 7 ) equations: 

dim T / Q \ 

53 £e s (X)a 9 (0 i ,^)U(X)=V(0 i ) .7 = 1,..., dim J 7 (4) 

i=i \ f / =i / 

In many cases, the space J 7 has a large dimension, so as to represent many functions of the 
continuous function space with a great precision, and the system (4) is very large (although one 
can generally choose T and (fa) so as to produce a very sparse system). 

When (3) has to be solved for many values of X (the so-called many query context), the reduced 
basis method can be used to speed up the overall computation, which is split into two parts. In 
the offline phase, we choose a linearly independent family B = {Ci, • • • ,Cn} of n vectors in J 7 , 
and compute and store the Q n-by-n matrices of each a q form [ie. the matrices A q whose 
coefficient is a q (Q, Q)) and the n-vector representing ip (ie. the vector whose zth coefficient is ip(0)) 
in the basis B (called the reduced basis). The offline phase does not depend on a particular value 
of X and can be done once. Then, for each value of X for which it(X) has to be computed, one 
can proceed to the online phase: using information stored during the offline phase, the following 
n-hy-n linear system is assembled and solved for (Mj(X))j = i i ... jn : 

n I Q \ 

53 53 e 9 (xK(c>, q «*(x) = ^(0) j = 1, • • • , n (5) 

i=l \q=l J 

n 

Then w(X) « m(X) is recovered by using w(X) = 53{tj(X)Cj- In many cases, {tt(X); X G T>} lies 

in a manifold of dimension much smaller than dim J 7 , and it is possible to choose a linear space of 
approximation of dimension n <C dim J 7 and thus, solve (5) much faster than (4) while keeping u 
sufficiently close to u. At the end of this section we will see a method to automatically choose an 
"effective" reduced basis, that allows accurate representation of m(X) for X G V. 
When the model output /(X) is a linear functional /(it(X)) of w(X), the surrogate output can be 
defined as: 

/(X):=/(n(X)) = 53n J (X)/(0) (6) 

9=1 
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whose parameter-independent reals f{d) (i = 1, . . . ,n) can be calculated and stored during the 
offline phase, allowing evaluation of /(X) without explicitly forming w(X), and leading to a meta- 
model whose complexity of evaluation depends only on its dimension n (and on Q) - and no more 
on the dimension of the original model dim T . 

1.2.2. Error bound 

An interesting feature of the reduced basis approach is that it comes with a provable error bound 
£ U (X) fully computable with a complexity independent of dim J 7 Nguyen et al. (2005). This error 
bound satisfies 

||u(X)-«(X)|| <e u (X) VXgD 

for a chosen Hilbert space norm ||-|| on T '. To present the error bound, we assume, for simplicity, 
that the a g 's are symmetric bilinear forms, and that a(-, -;X) is uniformly coercive, that is, a(X) 
defined by: 

a(X) = sup a(v, u;X) 

v£j r ,\\v\\=l 

satisfies a(X) > for all X G T>. 
We claim that: 

IKX)-«(X)||< MSIk vxgd 

where r(X) is the residual linear form, defined by: 

r(X)(u) = ip(v) - a(u(X),v; X) 



and W'W-p, is the dual norm on J 7 : 



\ F , = sup £(v) 

v &T, \\v\\=X 



Efficient procedures have been developed for efficient offline-online computation of ||r(X)||j-,, and 
a lower bound 5(X) < a(X), leading to a computable error bound on u: 



5(X) 

This error bound on u can be used to develop an error bound on the output. For example, if 
/(X) = /(u(X)) and /(X) is the surrogate output defined in (6), one can use 

e(X) = \\f\\ r e u (X) (7) 

which satisfies: 

/(X)-/(X)|<e(X) VXeD (8) 

as error bound on the output. 

1.2.3. POD-based procedure for reduced basis choice 

We now describe a way of selecting a reduced basis {Ci, ■ ■ ■ , Cn}- 

We randomly choose a finite subset H = {X 1 , . . . , X m } C P, and compute w(X) for each X G S. 
We put the coordinates of w(X) with respect to the basis {0i, . . . , 0dim t} of J 7 as columns of a 
snapshot matrix S: 



S 



( u{Tt% n(X 2 ) x ... M (X m ) x \ 
u(X% u(X 2 ) 2 ... u(X m ) 2 

U^toJ ^(X 2 ) dim ^ ... «(X m ) dim ^y 



(9) 



We now proceed with the Proper Orthogonal Decomposition (POD) of the S matrix: we compute 
{zi, . . . ,z n }, where Z{ is an eigenvector associated with the i th largest eigenvalue of the m-by-m 
symmetric matrix S T £IS (where 0, is the matrix of the scalar product <, > associated with ||-||, 
with respect to the {0i, . . . , 0dimj r } basis), and define the vectors of the reduced basis to be: 



One can show that the {£1, . . . , ( n } are solutions of the following optimization program: 

Minimize ^ ||w(X) — 7t[m(X)]|| 2 , under the constraints < &,Q >= 
xe~ 

where 7r is the orthogonal projector onto Span(£i, . . . , ( n ). 
Proper orthogonal decomposition (also known as Principal component analysis (PCA), or Singular 
value decomposition (SVD)) Chatterjee (2000), and variants of POD, are widely used in model 
reduction without error bounds Bui-Thanh et al. (2007); Bergmann and Iollo (2008). We also 
showed in Janon et al. (2010, submitted.) that POD reduced bases are efficient with respect to the 
obtained error bounds. 

1.3. Estimator on reduced model 

Using the reduced model to perform the sensitivity analysis is straightforward: replace every call 

to the full model / by a call to the reduced one /. This gives rise to an estimator Si converging, 
for N — > +oo, to the true value of the sensitivity index on /: 

~ _ VarE(F|X t ) 
VarF 

Sampling error of this estimator can be assessed using bootstrap Efron et al. (1993), see Algorithm 
1. 

However, doing so does not take in account the gap between Si and Si. If the metamodel is "too 
far" from the original model, the (1 — a)-confidence interval estimated using it will not contain the 
true value of Si with probability close to 1 — a. On the other hand, a moderate-fidelity metamodel 
might be well-suited to give a "rough" estimate of sensitivity indices - in some cases such a rough 
estimate is sufficient - but the user would like be informed that the metamodel he uses gives him 
a limited-precision estimator; such a limited precision would reflect in the increase in the width of 
the output confidence interval. 

2. Quantification of the two types of error in index estimation 

We now present our method for estimating the two types of error that occur in Monte-Carlo 
sensitivity index estimation on a reduced-basis metamodel. In the first part 2.1, we review the 
bootstrap, which we will use for the treatment of sampling error. In the second part 2.2, we show 
how to use reduced-basis bounds to assess metamodel error. 

2.1. Sampling error : bootstrap confidence intervals 

Sampling error, due to the Monte-Carlo evaluation of the variances in (1), can be quantified 
through an approximate confidence interval calculated using bootstrap Archer et al. (1997). 
We use the bias-corrected (BC) percentile method presented in Efron (1981); Efron and Tibshirani 
(1986). The principle of this method can be summed up the following way: let 9(X\, . . . , X n ) be an 



1 if i = j 
else. 
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estimator for an unknown parameter 9 in a reference population V . To get a point estimate of 6, 
one takes a random i.i.d. n-sample {x\, . . . , x n } from V, and computes 9(xi, . . . , x n ). In (nonpara- 
metric) bootstrap we repeatedly, for b = 1, . . . , B, sample . . . ,x n [b}} with replacement from 

the original sample {x\, . . . , x n } and get a replication of 9 by computing = 9(xi[b], . . . , £„,[&]). 
This way we get a sample TZ = {9[1], . . . , 9[B]} of replications of 9. 

Now see how this sample can be used to estimate a confidence interval for 9. We denote by $ the 
standard normal cdf: 

= — / exp 




and by <J> 1 its inverse. 

Using 1Z and the point estimate 9 = 9(x±, . . . , x n ), a "bias correction constant" zq can be estimated: 



z = $ 



-i 



'#{£[&] e ft s.t. §[b] < dy 



\ B J 

Then, for (3 g]0; 1[, we define the "corrected quantile estimate" 

q(/3) = <$>{2T + z p ) 

where zp satisfies Q(zp) = (3. 

The central BC bootstrap confidence interval of level 1 — a is then estimated by the interval whose 
endpoints are the q(a>/2) and q(l — a/2) quantiles of 1Z. 

This confidence interval has been justified in Efron (1981) when there exists an increasing trans- 
formation g, z G R and a > such that g{9) ~ N{9 — z a, a) and g{9*) ~ N{9 — z a, cr), where 
9* is the bootstrapped 9, for fixed sample {xi, . . . ,x n } and (hence) fixed 9 = 6(xi, . . . ,x n ). In 
practice, due to the complex analytic expressions (20) of the estimators we are going to bootstrap, 
it seems hard to prove that such a g exists. However, we give in Section 4.4 empirical evidence 
that, for the two estimators defined at (20), g can approximatively be chosen as identity. 
The full computation method of a BC bootstrap confidence interval for Si is given in Algorithm 1. 
The key advantage of bootstrapping our sensitivity estimators is that we do not require supple- 
mentary model evaluations to estimate a confidence interval; hence the computational overhead 
for getting a confidence interval (versus pointwise estimation only) remains quite modest. 

2.2. Metamodel error 

For a pair of samples ({X fc } fc=1 ^ {^-' k }k=i,...,N), we can use our metamodel output / and our 



metamodel error bound e to compute, for k = 1, . . . , N: 

y k = /(X fc ), y' k = f(X[ k , . . . , X*,, Xl X* lf ...,**) 

and: 

E k — E{1^ ), E k — E{Ji 1 , . . . , A^ 1? Aj , A i+1 , . . . , A p ) 

In this section, we find accurate, explicitly and efficiently computable bounds S™ and S^ 1 , de- 
pending only on yk,y' k ,Ek and e' k so that: 

S™ < Si < Sf (10) 

In other words, we want lower and upper bounds on the full model based sensitivity index estimator 

Si computable from surrogate model calls. 

Let: 

N 

R(a; y, fi, //') = ^{y' k - (a(y k - fi) + n')f 
k=i 
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where y = (y x , . . . , y N , y[, . . . , y' N ) and ll, //' G R 

By setting first derivative of R with respect to a to zero, making use of the convexity of R(-; y, y, y 1 
and using: 

jf J2k=i Vkv'k ~~ {jj ^2k=i Vkj yjf J2k=i y'k) 



S> 

one easily shows that: 



Ar J2k=l Dk \N ^k=l Vk 



Si = argmin R(a; y, y, y') 

aeR 



i N i N 

where: y = — ^ y k and y 7 = — ^ Z/fc- 

iV fc=L iV fe=l 

In other words, Si is the slope of the linear least squares regression of the {y' k }k on the {yk}- 
Define: 

Rmf(a\y,e,n,n') = ^2{ _ _ inf _ _ (z' - (a(z - //) + n')f \ (11) 

fe=1 ^G[j/fc-efe;yfc+£fc],z'G[^-e' fc ;^+4] J 

and: 

R sup (a;y,e,LL,Li') = ^ < _ sup _ (z 1 - (a(z - ll) + ju')) 2 > (12) 

fc=l [£e[y fc -£fe;yfc+e fc ],z'e[:y^-£' fc 'j4+ e 'fc] J 

where y = (yx, . . . , y N , y[, . . .,y' N ), e = (e u . . .,e N ,e[, . . -,e' N ). 
It is clear that: 

R inf (a;y,£,Li,Li') < R(a;y, n, [i') < R sup (a;y,s,Li,n') Va,//,//eR (13) 

Note that R, Ri n f and R sup are quadratic polynomials in a. We name a, (3, 7, ai n f, (3i n f, Jinf, o^ sup , Psu P 
and '-fsup their respective coefficients. In other words, we have: 

R(a; y, fx, //) = aa 2 + (3a + 7 

Rinf (a; f, £, fJ>, aO = <y inf a 2 + /3 in/ a + 7; n/ (14) 

R sup (a; y, e, At, //) = a SMp a 2 + (3 sup a + 7sup (15) 

These coefficients depend on fi and // . We do not explicitly write this dependence until the last 
part of our discussion. 

Using (13) we see that the quadratic function of a: 

(a inf - a)a 2 + (Jt inf - (3)a + j inf - 7 

is negative or zero; hence it takes a non-positive value for a = 0, and has a non-positive discrimi- 
nant: 

Tin/ " T < (16) 
(An/ - P? < 4(a m/ - a)( 7in/ - 7) (17) 

As (f3inf — (5) 2 > 0, Equations (16) and (17) above imply that ai n f — a < 0, and that: 

Pin} — Sinf < ft < An/ + ^in/ 



1 &s well on y (for a, ,8,7) and y and e (for the other coefficients) 
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for 6 in f = 2yJ(a inf - a)(7 in / - 7)- 

We now suppose that > 0. As a^f is computable from yk,y' k ,£k and e' k , one can practically 
check if this condition is met. If it is not the case, our bound can not be used. We expect that if 
the metamodel error is not too large, we have a in f « a and, as a > 0, the hypothesis a in f > is 
realistic. 

So, under this supplementary assumption, we have: 

in r 

argmm R(a; y, //, /x ) = -— > 

Now using the second part of (13) and the same reasoning on the non-positive quadratic function 
of a: R(a; y, fi, //) - R sup (a; y, e, //, //), we get that: a < a sup , and: /3 sup - 5 sup < /3 < (3 sup + 5 sup . 
Hence, 



argmin R(a; y, /i, //) < 



Psup 0~sup 

2ct sup 



where 5 sup = 2\J (a — a sup ){^ — j sup ). This comes without supplementary assumption, because 
o-sup > ol and a > 0, as the minimum of i?(-; y, /x, //) exists. 



As we clearly have 6i n f and 5 sup less than (or equal to) 5 := 2y(aj n / — a sup )(7i„/ — 7 sup ), we 
deduce that: 

Pinfiv, »') + 5(n, //') , /3 sup (/i, //) - 5(/x, //) 

\ argmin rtya, y, /1, /x j ^ 



2ttj n j(/x, 2a sup (/i, /x') 

where we have explicited the dependencies in /x and //. 
To finish, it is easy to see that we have: 

V:= [y-e;y + e] By (18) 

and: 

V' := [W-e';W + s] By' (19) 
(where y,y',e and e' denote, respectively, the means of (yk) k , iy'k)k > an d so that: 



mm_ 



< Si = argmin i*(a; y, //,//) < max 

2a inf {fj,, fi') J a nev^'eP \ 2a sup ([i, /x'J 

Hence, (10) is verified with: 

S™= min L /M/^O + WA § „ = ^ ( j sup ^^')-5^^) \ 
nep,ii'ep' \ 2a in f([x,n') J pepper' \ 2a sup (fi, /x'J J 

It is clear that S™ and Sf 1 are computable without knowing the y k s and y' k s. 
In practice, we compute approximate values of S™ and Sf 1 by replacing the min and max over 
V x V by the min and max over a finite sample ScPx V. See Algorithm 2 for a summary of 
the entire computation procedure. 



3. Combined confidence intervals and parameters choice 

3.1. Combined confidence intervals 

In the last section, we have seen how to separately assess sampling error and metamodel error. To 
take both error into account simultaneously, we propose using bootstrap confidence intervals (see 



9 



Section 2.1) by calculating B bootstrap replications of S™ and S^ , where, for b — 1, . . . , B each 
bootstrap pair (s?[b];SV[b]) is computed using {yk)k£L b ' (y'k)keL b as surrogate output samples, 
and associated error bounds (sk)keL > (^fc)fcei, > wnere L b is a list of N integers sampled with 
replacement from {1, . . . , N}. 

The BC bootstrap confidence interval procedure (see 2.1) can then be used to produce a 1 — a-level 
confidence interval [S™ a , 2 ; S^{_ a i 2 ) f° r ^Ti an d a confidence interval [S^ a / 2 ; S^[_ a , 2 ] for Sf 1 . We 
then take as combined confidence interval of level 1 — a for Si the range [S™ a / 2 ; S^i-a^]- This 
interval accounts for sampling and met amo dels error simultaneously. 

e sampling 

Optionally, we can introduce a postulated uncertainty on the error bounds through what we 
call e sampling. In e sampling, the 6 th bootstrap replicates for S™ and Sf 1 are computed using 
( £ t)k£L ' \ £ '*k) as error bounds, where e* k and e'* k are sampled independently from a uniform 
distribution in [rjkEk'i £ k] an d [Vk £ k'i e 'k\ > where r] k ,r]' k e [0; 1] are alleged effectivities of our error 
bound, that is, an indicator of the ratio between the true errors \yk — Vk\ an d Sk (and between 
Wk ~ y'k\ an d £fc)- Setting effectiveness close to zero narrows confidence intervals, putting more 
trust in the reduced model than in the error bound, which is considered too pessimistic; on the 
contrary, effectiveness close to one means that error bound does not overestimate true error too 
much and that the error can not be considered too smaller than it. 
The procedure for obtaining confidence intervals is summed up in Algorithm 3. 

3.2. Choice of reduced basis size and Monte- Carlo sample size 

When doing a Monte Carlo estimation of sensitivity indices using a reduced basis metamodel, by 
means of confidence intervals computed with the strategy described above, one has to choose two 
important parameters : the sample size (N) and the number of elements in the reduced basis 
(n). Increasing N and/or n will increase the overall time for computation (because of a larger 
number of surrogate simulations to perform if N is increased, or, if n is increased, each surrogate 
simulation taking more time to complete due to a larger linear system to solve). However, increase 
in these parameters will also improve the precision of the calculation (thanks to reduction in 
sampling error for increased N, or reduction in metamodel error for increased n). In practice, one 
wants to estimate sensitivity indices with a given precision (ie. to produce (1 — a)-level confidence 
intervals with prescribed length), and has no a priori indication on how to choose N and n to do 
so. Moreover, for one given precision, there may be multiple choices of suitable couples (N,n), 
balancing between sampling and metamodel error. We wish to choose the best, that is, the one 
who gives the smallest computation time. 

The aim of this section is to describe a simple computational model that helps us in making a good 
choice of sample size and reduced basis size to produce a confidence interval of a desired precision. 

Formulation as a constrained optimization problem 

On the one hand, we evaluate computation time: an analysis of the reduced basis method shows 
that the most costly operation made during an online evaluation (see Section 1.2) is the resolution 
of a linear system of n equations; this resolution can be done (e.g., by using Gauss' algorithm) 
with 0(n 3 ) operations. This has to be multiplied by the required number of online evaluations, 
i.e. the sample size N. Hence, we may assume that computation time is proportional to N x n 3 . 
On the other hand, the mean length of the (1 — a)-level confidence intervals for Si, . . . , S p can be 
written as the sum of two terms. The first, depending on N, accounts for sampling error and can 
be modelled as 




10 



for a constant Z a > 0. The assumption of 1/v N decay is heuristically deduced from central limit 
theorem. 

The second term, which accounts for metamodel error, is assumed to be of exponential decay 
when n increases: C/a n , where C > and a > 1 are constants. This assumption is backed up by 
numerical experiments as well as theoretical works Buffa et al. (2009). 

Once this analysis has been done, we translate our problem into the following constrained mini- 
mization program: 

Find (N*,n*) = argmin n 3 x iV so that + — = P (21) 

(n,jV)eR + xR + a n 

where P is the desired average precision for the confidence intervals. 

Note that we converted the discrete design variables N and n to continuous positive variables so 
as to use the standard tools of continuous optimization; once optimum of the continuous problem 
have been found, we just round it to the nearest integer couple to recover a near-optimal integer 
solution. 

Resolution of the optimization problem 

The constraint in (21) is equivalent to the conjunction of the following two equations: 

N -(-^s) (22) 

P ~ ~ n > ( 23 ) 



a 

so that the function to minimize over I = ]n c ; +oo[ (where n c = ln(C/P)/ln(a) has been chosed 
so as to satisfy (23)) is: 

(j){n) = (Z a ) —2 

This function is differentiable on I and tends to +oo as n — > n c and n — > +oo; hence it has a 
minimizer n* G / that satisfies <j)'(n*) = 0, which is equivalent to: 

3 (24) 



On J, (24) is equivalent to: 



Pa n * -C 2C In a 
3P 



Ti — — a n 



2C\na 2C\na 



3P 

Now let ib be the function defined on R by ib(x) = x — — a x . By remarking that ib'(x) 

26 ma 



1 — - — a x is continuous and nonzero on 
2C 



2(7 \ 

In ( — — ) / In a; +oo 



D I, one has that ib is injective on 



3PJ 

I, and so (24) has at most one solution in /. Thus (24) has exactly one solution in /, of which 
an approximate value can be found by using bisection method Press et al. (1992) on [n c + e; L\ 
where e > is small enough and L is a sufficiently large. Once n* has been found, we can find the 
optimal N* by setting n = n* in (22). 
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Estimation of the parameters 

The last question that remains to address is the estimation of the Z a , a and C constants. The Z a 
parameter is estimated by running the estimation procedure on the metamodel, for fixed N and 
n, estimating combined BC bootstrap confidence intervals (Section 2.1) and taking: 

Z a = ^/N {pifl-a/2 ~ Si* + ^>/2 — 

where the factor multiplying Vn is the estimated "Monte-Carlo part" of the error. 
The a and C parameters are estimated by running an estimation procedure, for a single N fixed, 
and different reduced basis sizes rii, . . . , tik, and measuring, for each reduced basis size, the average 
metamodel error e(rifc) = — S™, for k — 1, . . . , K. The { e(rii)); . . . ; (n^, e(%))} pairs are 
then used to fit the exponential regression model e(n) = C/a n . 

If one wants to estimate the sensitivity indices with respect to all variables i — 1, . . . ,p for a single 
value of N and n, one can use bootstrap procedure to estimate Monte-Carlo errors E±, . . . , E p for 
each of the p sensitivity indices estimators: 

— °i,l-a/2 ~ D i i* °i,a/2 ~ D i 

and then to take: 

Z a = >J 

4. Numerical results and discussion 

In this section, we test our combined confidence interval procedure described earlier, and compare it 
with Monte-Carlo estimation on the full model (with bootstrap to assess sampling error), and with 
the procedure described in Storlie et al. (2009) and implemented in the CompModSA R package. 
Our criteria of comparison are the CPU time needed to compute the intervals and the lengths of 
these intervals (the smaller the better). 

In all our tests we take a = .05 and B = 2000 bootstrap replications. 
4-1. Model set-up 

Let u, a function of space x G [0; 1] (note that space variable x is unrelated to input parameter 
vector X) and time t e [0, T] (T > is a fixed (i.e., known) parameter) satisfying the viscous 
Burgers' equation: 

du 19.,. d 2 u . 

Tt + 2 *<" > - "a? = * (25) 



where v G Rj~ is the viscosity, and ip G C° (JO, T], L 2 (}0, l[)j is the source term. 
For it to be well-defined, we also prescribe initial value Uq G -ff 1 (]0, 1[): 

u(t = 0, x) = u (x) Vx G [0; 1] (26) 

and boundary values bo,bi G C ([0,T]): 



u(t,x = 0) = b (t) 
u(t,x= 1) = h(t) 



Wt G [0; T] (27) 



Where 6o, 6i and uq are given functions, supposed to satisfy compatibility conditions: 

u (0) = b (0) and u (l) = 6^0) (28) 
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The initial u$ and boundary values bo and b\ are parametrized the following way: 



n(bo) n(bi) 

bo(t) = b 0m + K° sin( W f t) b x {t) = b lm + ]T A b ^ sin(^t) 

1=1 1=1 

nr(/)ns(/) n(ito) 

= f m + J2 E <sin(w/ T t) sin(wf x) u Q (x) = (u 0m ) 2 + £ sin(^°x) 

J=l p=l Z=i 

The values of the angular frequencies wf°, wf 1 , w/ T , Wp 5 and u^°, as well as their cardinalities 
n(bo), n(pi), riT(f), ns(f) and n(uo) are fixed (known), while our uncertain parameters, namely: 
viscosity u, coefficients b 0m , b lm , f m and u 0m , and amplitudes (A^ ) f^f 1 ), , o 

\ / (=l,...,n(oo) \ / (=l,...,n(bi) 

(Af„) and , „/. s live in some Cartesian product of intervals V' , subset 

V l PJl=l,...,n T (f);p=l,...,n s (f) K 1 n=l,...,n{u ) 

Of R, 1+4+n ( b °) +n ( 6l )+™T(/)ns(/)+«(«o) 

However, the compatibility condition (28) constraints &om and 6i m as functions of the other pa- 
rameters: 

n(uo) 

b 0m = (u 0m ) 2 and b lm = (u 0m ) 2 + £ A? sin(^ u °) (29) 

1=1 

so that the "compliant" uncertain parameters actually belong to V defined by: 



— |X — (V, &0m, ^1°, • • • , ^n°(feo)' ^ lm ' ^l 1 ' ' ' ' ' ^n(6i)J ^11 > ^12, • • • j ^l,n s (/)> 

■ ■ ■ , 4n s( /)' • • • , A Lu),nsUV M °™> ' • ' ' A n(n )) e p/ satisfying (29) \ (30) 



In Janon et al. (2010, submitted.) , we gave an example with many more parameters. To illustrate 
our sensitivity analysis methodology without overloading the text, we choose an example with a 
reduced number of parameters. 

The solution u = w(X) depends on the parameter vector X above. 

The "full" model is obtained by discretizing the initial-boundary value problem (25), (26), (27), 
using a discrete time grid {tk = kAt}k=o,...,T/At, where At > is the time step, and, space-wise, 
using P 1 Lagrange finite elements built upon an uniform subdivision of [0; 1]: {xi = i/ftf}, for 
i = 0,... ,7V. Our full output is: 

/(X) = f(u(X)) = i £ u(t = T, x = Xl ) 
■ /v i=0 

The reduced basis method is then applied to yield a surrogate solution u of (25), (26), (27), as 
well as an error bound e u on u. Due to non-linearity and time-dependence of (25), as well as 
parametrization of the boundary values, the reduced basis methodology is not as simple as the one 
presented in Section 1.2 of the present paper. The reader can refer to Janon et al. (2010, submitted.) 
for full details on discretization and reduction of this model. Error bound on output e is obtained 
by following (7). 

In our numerical experiments, we take M = 60, At = .01, T = .05, ns(f) = n?(f) = n(bo) = 

n(6i) = 0, n(u ) = 1, up = 0.5, A\° = 5 and f m = 1. 

Reduced basis are found using POD-based procedure with #H = 30. 

The two input parameters are assumed independent and uniformly distributed. The table below 
contains the ranges for them, and also the "true" values of the sensitivity indices, which have been 
calculated (in more than 14h CPU time) using a Monte-Carlo simulation with large sample size 
JV = 4x 10 6 (so as to BC bootstrap confidence intervals of length < 10 -2 ) on the full model: 



13 



nu 



Figure 1: Output / of the Burgers' model, plotted as a function of v and mq 



Parameter 


Range 


95% confidence interval for sensitivity index 


V 


[1 ; 20] 


[0.0815;0.0832] 


UQm 


[-0.3 ; 0.3] 


[0.9175;0.9182] 



The output, as a function of the two uncertain parameters v and u 0m is plotted at Figure 1; as 
one can see it is nonlinear with respect to the input parameters. 

4-2. Convergence benchmark 

Figure 2 shows the lower S m and upper S M bound (defined in Section 2.2) for different reduced 
basis sizes (hence different metamodel precision) but fixed sample of size N = 300, as well as 
the bootstrap confidence intervals computed using the procedure presented in Section 3.1. This 
figure exhibits the fast convergence of our bounds to the true value of the sensitivity index as the 
reduced basis size increases. We also see that the part of the error due to sampling (gaps between 
confidence interval upper bound and upper bound, and between confidence interval lower bound 
and lower bound) remains constant, as sample size stays the same. 

4 ■ 3. Choice of n and N 

We now discuss the numerical results obtained when using the parameter tuning procedure (Section 
3.2). We have done "pre-runs" for N = 300, and different reduced basis sizes {7, 8, ... , 12} of the 
combined confidence interval procedure, to yield the following estimates: 

C = 197.69 a = 2.789 Z 05 = 2.6407 

To assess validity of this estimation, and to check our modelisation of the bound precision and the 
execution time, we plot the cube root of the CPU time (Figure (3)) and the precision of the bound 
defined in Section 2.2 (Figure (4)). We can check that these hypotheses are reasonably satisfied. 
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0.4 



0.35 



0.3 



0.25 



0.2 



0.15 



0.05 



-0.05 



1.05 



lower bound S™ hat 
upper bound S, hat 
lower bound, 95% CI S, 05/2 ™ hat 
upper bound, 95% CI S, ^ ' 05/2 hat 



9 10 
reduced basis size 
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lower bound S™ hat — 
upper bound S; hat -- 

lower bound, 95% CI S, 05/2 ™ hat 

upper bound, 95% CI S, ^ ' 05/2 hat 















0.75 



0.95 



0.9 



0.85 
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Figure 2: Convergence benchmark for sensitivity indices estimation in the Burgers' model, top: variable v, bottom: 
variable uo m - We plotted, for a fixed sample size TV = 300, estimator bounds S m and S M defined in (2.2), and 



endpoints «S™ 02 5 and S, 



of the 95% confidence interval, for different reduced basis sizes. 
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0.92 




0.72 1 1 1 1 1 1 

7 8 9 10 11 12 

reduced basis size 



Figure 3: Cube root of the CPU time necessary to do estimations, as a function of the reduced basis size of n. 
Section 3.2 assumes this function is linear. 




7 8 9 10 11 12 

reduced basis size 

Figure 4: Line: reduced basis "precision", i.e. mean Sf 1 — S™, as a function of reduced basis size; dashes: result of 
the fit of an exponential regression model: C/a n . 
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Precision p 


Reduced basis size n* 


Sample size N* 


0.005 


12.4437 


354491 


0.02 


11.1095 


22057.6 


0.05 


10.0501 


3698.95 


0.08 


9.59689 


1442.7 


0.09 


9.48332 


1139.47 



Figure 5: Optimal reduced basis and sample sizes calculated using the strategy described in Section 3.2. 




Figure 6: Normal empirical quantile-quantile plots of the distributions of S™ (left) and S v (right). 



One can find in Figure 5 the computed optimal reduced basis sizes n* and sample sizes N* using 
resolution of the optimization problem (21), for various precisions p. 

All these values have been computed in 5.77 s CPU time, including the time necessary to estimate 
C, a and a. 

To check for the efficiency of this parameter tuning strategy, we choose a target precision of p = .02. 
In the table 5, we read that we should take n ~ 11 and iV « 22000. 

Conducting the combined confidence interval estimation with these parameters give intervals 
[0.0659997; 0.0937285] for sensitivity index for u, and [0.914266; 0.926452] for sensitivity with re- 
spect to Mom- These confidence intervals have mean length: 

- (0.0937285 - 0.0659997 + 0.926452 - 0.914266) = 0.0199575 « 0.02 

as desired. 

This computation took 52 s of CPU time to complete (including a metamodel offline phase of 1 s). 
4-4- Normality of the bootstrap distributions 

We give in Figure 6 the empirical normal quantile-quantile plots of the bootstrap replicates 
{S?[l],...,S?[B]}™d 

As these plots are close to a line, the bootstrap distributions are approximately normal. 
4-5. Optimality of our metamodel error bound 

We checked for near optimality of the metamodel error bound 2.2 by comparing it with the values 
of the optimization problems: nun ip and max ip where: 

, . jf ^k=i UkUk ~ [jf SfcLi Vk) {jj SfcLi y'k) 
' y = £EiL(ifc) 2 -(*£iLy fc ) 2 
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Surface response (metamodel) 


Mean confidence interval length 




CPU time 


qreg: quadratic regression 


0.081 


0.996799 


143.59 s 


mars: multivariate adaptive regression splines 


0.075 


0.9998506 


218.716 s 


our approach 


0.019 


N/A 


52 s 



Figure 7: Results of CompModSA's sensitivity function on our model, for two fitted response surfaces. R 2 is an 
indicator of the metamodel fit (values close to unity suggest good fit). The last line recalls the results of the 
experiment in Section 4.3. 



and: 

N N 

y = II iVk - £k, Vk + £fc] x J] [y' k - 4; y' k + e' k ] 

k=l k=l 

These problems, of large dimension 27V, give the optimal values of S™ and Sf 4 satisfying (10). They 
can be solved with simulated annealing Kirkpatrick et al. (1983); Pardalos and Romeijn (2002). 
Our bound gave results very close to the optimal ones, for a smaller computational cost than using 
simulated annealing. 

4-6. Comparison with estimation on the full model 

To obtain a result of the same precision, we carry a simulation for N = 21000 (sample size can 
be chosen smaller than before, as there will be no metamodel error) on the full model; we get the 
bootstrap confidence interval with mean length of ~ 0.0193. 

This computation takes 308 s of CPU time to complete. Hence, on this example, using a reduced- 
basis surrogate model roughly divides overall computation time by a factor of 5.9, without any 
sacrifice on the precision and the rigorousness (as our metamodel error quantification procedure 
is fully proven and certified) of the confidence interval. We expect higher time savings with more 
complex (for example, two- or three-dimensional in space) models. 

4-7. Comparison with CompMoalSA 

We compared our results with the ones obtained using the R CompModSA version 1.2 package 
downloaded at compmodsa. This package implements the method described in Storlie et al. (2009) 
for assessing sampling error as well as metamodel error. It does not make use of the reduced basis 
output error bound, but uses a non-intrusive method to fit a metamodel using a reasonable number 
of full model evaluations. 

We fed into CompModSA procedure 50 such full model outputs (which took 0.22 s CPU to com- 
pute). We then tried various non-intrusive metamodels {surface responses), and reported the 
results into Figure 7. We used as parameters: n.mc.T=0 (we do not want any total index com- 
putation), n.mc.S=23000 (sample size), n.samples=l (one run), n.CI=300 (generate confidence 
intervals using 300 bootstrap replications). We contributed 2 the option CI.S, which is set to 
TRUE to compute bootstrap confidence intervals for the main effect index. 

By looking at results in Figure 5, we can see that in this case, our approach is clearly superior, 
both in terms of precision and computation time. To achieve this result, we took advantage of 
the particular formulation of the original model which allows the reduced basis methodology to 
be efficiently applied; CompModSA, due to its non-intrusive nature, is easier to use on a generic 
"black box" model. 



2 patch available at http://ljk.imag.fr/membres/Alexandre.Janon/compmodsa.php 
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Conclusion and perspectives 



We presented a methodology to make a rigorous quantification of the impact of the sampling error 
and the metamodel error on the sensitivity indices computation, when a reduced-basis metamodel 
is used. Sampling error is handled by a classic bootstrap procedure, while metamodel error is 
managed using a new bound on the sensitivity index estimator. Quantification of those two types 
of errors permits not only a certification on the performed estimation, but also gives a way to tune 
the optimal parameters (reduced basis and optimal sample sizes), for a given desired precision on 
the indices. We have shown on a concrete example the superiority of this method when compared to 
the use of the full model, or non-intrusive (quadratic regression, MARS) metamodels. Our method 
can be applied to other Monte Carlo (or quasi Monte Carlo) estimators, and to other metamodels 
which provide an error bound similar to the one provided by the reduced basis framework. 

Acknowledgements. We wish to thank Jean-Claude Fort for a suggestion which we exploited to 
perform our computation of the metamodel-induced error. We also thank Anestis Antoniadis 
and Ingrid Van Keilegom for advice on bootstrap methodology. - This work has been partially 
supported by the French National Research Agency (ANR) through COSINUS program (project 
COSTA-BRAVA n°ANR-09-COSI-015). 

A. Algorithms 



Algorithm 1: 

1. Draw from X distribution two independent samples of size N: {X fc } and {X' fc }. 

2. Tabulate necessary model evaluations: for k — 1, . . . , N: 

(a) set X «- X fc ; 

(b) compute y k = /(X); 

(c) swap Xi and Xf] 

(d) compute y' k = /(X); 

(e) swap Xi and Xf back. 

3. Compute Sf 

= '/ \~2 

N 2^k=lilk \n 2^k=iykJ 

4. Repeat, for b = 1, . . . , B: 

(a) Draw at random a list L of length N, with replacement from {1, . . . , N}. 

(b) Compute replication Si[b] : 

- ril n J2keLVky'k — (j;12k£Lyk) (jfJ2keLy' k ) 

S i[°\ = ,2 

TV SfeeL Vk ~ \n SfceL Uk ) 

5. Compute z : 

Tq = $ _a f #{be{l,...,B} s .t. Sj[b] < Sj} ^ 
1 p / t 2 ' 



where $(z) = , — / exp ( | dt 
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6. Look up for z a /2 so that: 






®{z a/2 ) = a/2 




and take Zi_ a / 2 = —z a / 2 , 


satisfying: $(z 1 _ a / 2 ) — 1 — a /2- 




7. Compute q(a/2) and q(l 


-a/2): 




q(a/2) = 


$(2S^ + 2„/ 2 ), o(l — a/2) = $(2^n + a/2) 




8. Compute S i>a / 2 and S iy i_ 


^2, the q(a/2) and g(l — a/2) quantiles of {^[1], . 


;Si[B]}. 


9. Output Si jCt /2\ as confidence interval for Si of level 1 — a. 





Algorithm 2: 

1. Draw from X distribution two independent samples of size N: {X fc } and {X' fc }. 

2. Tabulate necessary model evaluations: for k = 1, . . . , N: 

(a) set X «- X fc ; 

(b) compute = /(X) and e k = e(X) ; 

(c) swap Xj and A"- fc ; 

(d) compute = /(X) and e' k = e(X); 

(e) swap Xj and X- fc back. 

3. Compute y, y', e and e', the respective means of {y~k}, {y'k}j { £ k} and {s' k }. 

4. Choose a finite subset H of the set V x V where P and P' are defined by (18) and (19). 

5. Repeat, for (ji, ji') G S: 

(a) By using (11) and (12), compute R in f(a;y,e, /i, //) and R sup (a; y, e, /i, //) for three 
different values of a; 

(b) deduce a inf , (3^,^, a sup , (3 sup and 7 sup satisfying (14) and (15); 

(c) if aj re y < 0, exit with failure; 



(d) compute 5 = 2^ (a inf - a sup ){^ in} - -f sup ); 

(e) compute: 

6. Output: 

= , nun AiO £f = f max (/i, p') 



Algorithm 3: 




1. Follow steps 1. and 2. of Algorithm 2. 




2. Compute bounds S™ and ^ M using steps 3.-6. of Algorithm 2. 




3. Repeat, for b = 1, . . . , B: 




(a) Draw at random a list L of length N, with replacement from {1, 


...,7V}. 


(b) If using e-sampling: for k G L, sample e* k uniformly in [rjkSk] £fc] 


and e'* k uniformly 


in Mb4; 4] • 
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(c) Else: take, for k G L, e£ — e k and e'* k = e' k . 

(d) Compute bounds S™[b] and S™[b] using steps 3.-6. of Algorithm 2, with {yk)k^z, 
instead of {Vk) k=ly .. tN , (y'k)keL instead of {y' k )k=i,...,N as s am ple data, and (e* k ) kl - L 
instead of (£fc) fc=1 N , and \^'*k) k L instead of (e' k ) k=1 N as error bounds. 

4. Compute, for w G {m, M}, the two bias correction constants: 



where <&(z) = — / exp dt. 

\/2lT J-oo \ 2 1 



t 2 ' 



B 



5. Look up for z a /2 so that: 

$(z a/2 ) = a/2 

and take Zi_ a /2 = —z a /2, satisfying: $(z 1 _ Q ,/ 2 ) = 1 — a/2. 

6. Compute § m (a/2) and q M (l - a/2): 

q m (a/2) = §(2zf + z a/2 ), q M (1 - a/2) = $(2^P + *l-«/ 2 ) 

7. Compute S™ /2 and S$_ a/2 , the g m (a/2) and g A/ (l-a/2) quantiles of {S^[l], S™[B}} 
and {Sf [1], . . . , Sf[B]}, respectively. 

8. Output S™ a i 2 ] Sf'[_ a i 2 as combined confidence interval for Si of level 1 — a. 



References 

Archer, C, Saltelli, A., Sobol, I., 1997. Sensitivity measures, ANOVA-like techniques and the use 
of bootstrap. Journal of Statistical Computation and Simulation 58, 99-120. 

Bergmann, J., Iollo, A., 2008. Numerical methods for low-order modeling of fluid flows based on 
POD . 

Boyaval, S., Bris, C, Maday, Y., Nguyen, N., Patera, A., 2009. A reduced basis approach for 
variational problems with stochastic parameters: Application to heat conduction with variable 
robin coefficient. Computer Methods in Applied Mechanics and Engineering 198, 3187-3206. 

Buffa, A., Maday, Y., Patera, A., Prud'homme, C, C, T., 2009. A priori convergence of the 
greedy algorithm for the parametrized reduced basis. Mathematical Modelling and Numerical 
Analysis . 

Bui-Thanh, T., Willcox, K., Ghattas, O., van Bloemen Waanders, B., 2007. Goal-oriented, model- 
constrained optimization for reduction of large-scale systems. Journal of Computational Physics 
224, 880-896. 

Chatterjee, A., 2000. An introduction to the proper orthogonal decomposition. Current Science 
78, 808-817. 

compmodsa, . CompModSA. http://www.stat.unm.edu/~storlie/CompModSA/. 



21 



Efron, B., 1981. Nonparametric standard errors and confidence intervals. Canadian Journal of 
Statistics 9, 139-158. 

Efron, B., Tibshirani, R., 1986. Bootstrap methods for standard errors, confidence intervals, and 
other measures of statistical accuracy. Statistical science 1, 54-75. 

Efron, B., Tibshirani, R., Tibshirani, R., 1993. An introduction to the bootstrap. Chapman & 
Hall/CRC. 

Grepl, M., Maday, Y., Nguyen, N., Patera, A., 2007. Efficient reduced-basis treatment of nonaffine 
and nonlinear partial differential equations. Mathematical Modelling and Numerical Analysis 
41, 575-605. 

Grepl, M., Patera, A., 2005. A posteriori error bounds for reduced-basis approximations of 
parametrized parabolic partial differential equations. Mathematical Modelling and Numerical 
Analysis 39, 157-181. 

Helton, J., Johnson, J., Sallaberry, C, Storlie, C, 2006. Survey of sampling-based methods for 
uncertainty and sensitivity analysis. Reliability Engineering & System Safety 91, 1175-1209. 

Homma, T., Saltelli, A., 1996. Importance measures in global sensitivity analysis of nonlinear 
models. Reliability Engineering & System Safety 52, 1-17. 

Janon, A., Nodet, M., Prieur, C, 2010, submitted. Certified reduced-basis solutions of vis- 
cous Burgers equations parametrized by initial and boundary values. Preprint available at 
http : //hal . inria . f r/inria-00524727/en. 

Kirkpatrick, S., Gelatt Jr, C, Vecchi, M., 1983. Optimization by simulated annealing. Science 
220, 671. 

Marrel, A., Iooss, B., Laurent, B., Roustant, O., 2009. Calculations of sobol indices for the gaussian 
process metamodel. Reliability Engineering & System Safety 94, 742-751. 

Nguyen, N., Veroy, K., Patera, A., 2005. Certified real-time solution of parametrized partial 
differential equations. Handbook of Materials Modeling , 1523-1558. 

Pardalos, P., Romeijn, H., 2002. Handbook of global optimization. Volume 2. Kluwer. 

Press, W., Flannery, B., Teukolsky, S., Vetterling, W., 1992. Numerical recipes in C: the art of 
scientific programming. Cambridge U. Press, Cambridge, England . 

Saltelli, A., 2002. Making best use of model evaluations to compute sensitivity indices. Computer 
Physics Communications 145, 280-297. 

Saltelli, A., Chan, K., Scott, E., . Sensitivity analysis. 2000. 

Sobol, I., 2001. Global sensitivity indices for nonlinear mathematical models and their Monte 
Carlo estimates. Mathematics and Computers in Simulation 55, 271-280. 

Storlie, C, Swiler, L., Helton, J., Sallaberry, C, 2009. Implementation and evaluation of non- 
parametric regression procedures for sensitivity analysis of computationally demanding models. 
Reliability Engineering & System Safety 94, 1735-1763. 

Veroy, K., Patera, A., 2005. Certified real-time solution of the parametrized steady incompressible 
Navier-Stokes equations: Rigorous reduced-basis a posteriori error bounds. International Journal 
for Numerical Methods in Fluids 47, 773-788. 



22 



